Digital Twins and Reproducibility

Bridging Point-and-Click COPASI Software with Open-Source and Flexible R Programming language

First Published: Mar 14, 2025

Last updated: mai 19, 2025


Bastien CHASSAGNOL |
Post-Doc
LPSM and Laval

Lilija Wehling |
Ph.D.
VPE BioMed X Institute

Atreyee Banerjee |
Ph.D.
Mainz Universität

Soria Gasparini |
Ph.D.
Heidelberg University

Sanjana Balaji Kuttae |
Computational Biologist (Research Assistant)
VPE BioMed X Institute


Slides available on GitHub

Talks’ Layout

  • Reproducibility Crisis. 👋
  • IBD ODE Model 🕘
  • IBD Numerical Simulations 🕘
  • BioModels database ⬇️

Reproducibility Crisis

Is there a Reproducibility crisis?


Classes of cell-cell communications cover:

Nature’s survey of 1,576 researchers on reproducibility in research. Reproduced from Fig1, 1,500 scientists lift the lid on reproducibility, from Baker (2016).

Nature’s survey of 1,576 researchers on reproducibility in research. Reproduced from Fig1, 1,500 scientists lift the lid on reproducibility, from Baker (2016).

Physicists and chemists show the most confidence, in contrast with social or biological studies. Reproduced from Fig2, 1,500 scientists lift the lid on reproducibility, from Baker (2016).

Physicists and chemists show the most confidence, in contrast with social or biological studies. Reproduced from Fig2, 1,500 scientists lift the lid on reproducibility, from Baker (2016).

Half of public mathematical models stored in BioModels are not curated. Results drawn from 17.05.2025, from Malik-Sheriff et al. (2020).

Half of public mathematical models stored in BioModels are not curated. Results drawn from 17.05.2025, from Malik-Sheriff et al. (2020).
  • \(52\%\) considering there’s a significant crisis

  • \(90\%\) considering there’s a Reproducibility crisis generally speaking

  • From Open Science Collaboration (2015) analyses, \(40\%\) of studies in psychology are trustworthy…

  • … and only \(10\%\) in cancer biology, from Begley and Ellis (2012)

Different types of reproducibility

Aka Methodological Reproducibility:

  • Ability to repeat the same experiment using the same methodology and obtain the same results
  • Example in scRNA-seq: 2 teams working on the same tissue, get similar distribution of cell type abundances
  • Reproduce statistical results derived from data analysis (\(p\)-values, effect sizes, confidence intervals, …) using the same statistical methods (or slightly similar ones)

  • Also named Inferential Reproducibility

  • Examples in sensitivity analyses: does a change of the random seed or the sampling process induce significant changes of biological findings?1

  • Also named Result Reproducibility, can encompass all previously mentioned reproducibility approaches.

  • Ability to reproduce the computational aspects of a study, which includes reproducing the figures, tables, or other outputs from the data and code provided.

  • Example in mathematical models: given the same fixed parameters, can you recover the same steady-state conditions?

Reproduce figures and main results from biological papers: Impossible mission?

Main factors contributing to irreproducible research. Unavalaible code and/or datasets is one of the most quoted reasons, while being one of the easiest to address with. Reproduced from Fig.4, Baker (2016).

Main factors contributing to irreproducible research. Unavalaible code and/or datasets is one of the most quoted reasons, while being one of the easiest to address with. Reproduced from Fig.4, Baker (2016).

Mentoring and training emerge out as the main incentives, along with more stringent journal guidelines. Reproduced from Fig.5, Baker (2016).

Mentoring and training emerge out as the main incentives, along with more stringent journal guidelines. Reproduced from Fig.5, Baker (2016).

Repro-hackathons as community-driven events to overcome reproducibility crisis?

Comparison of differentially expressed genes among 12 groups, illustrated by upset plots. Reproduced from Fig.6, Cokelaer, Cohen-Boulakia, and Lemoine (2023)

Comparison of differentially expressed genes among 12 groups, illustrated by upset plots. Reproduced from Fig.6, Cokelaer, Cohen-Boulakia, and Lemoine (2023)
  • Core objective: 12 group students aim at reproducing the main findings from a previously published study, namely the subset of DEGs.

  • Two-stage Hackathon:

    1. Free choice of differential analysis model, genome reference and filtering thresholds (scenario A)
    2. Constrained choice of methods and hyperparameters (scenario B)
  • Conclusion: majority of the differentially expressed genes were common to only a few groups (and no DEG found across all the 12 groups), while all pipelines were runnable.

Repro-hackathon to curate biological models to unified SBML standards

BioModels x BioQuant x BioMedX Team VPE sponsors.

BioModels x BioQuant x BioMedX Team VPE sponsors.
  • Objective: curate 11 mathematical models of human disease. Cross-disciplinary teams: one biologist, one data scientist, and one computational modeller.
  • Methods: reproduce at least one output (or recode when only hard-coded ODE equations were provided) of a given biological model + curate species to standard ontologies.
  • Results: 2 models were completely curated, 4 models were almost recoded and reproduced, 5 models were partially reproduced.

An ODE of the Immune System

Inflammatory Bowel Disease: an ODE model

IBD ODE network is depicted in Figure 1.

Code
flowchart LR
      A["$$x^2$$"] -->|"$$\sqrt{x+3}$$"| B("$$\frac{1}{2}$$")
      A -->|"$$\overbrace{a+b+c}^{\text{note}}$$"| C("$$\pi r^2$$")
      B --> D("$$x = \begin{cases} a &\text{if } b \\ c &\text{if } d \end{cases}$$")
      C --> E("$$x(t)=c_1\begin{bmatrix}-\cos{t}+\sin{t}\\ 2\cos{t} \end{bmatrix}e^{2t}$$")

 flowchart LR
      A["$$x^2$$"] -->|"$$\sqrt{x+3}$$"| B("$$\frac{1}{2}$$")
      A -->|"$$\overbrace{a+b+c}^{\text{note}}$$"| C("$$\pi r^2$$")
      B --> D("$$x = \begin{cases} a &\text{if } b \\ c &\text{if } d \end{cases}$$")
      C --> E("$$x(t)=c_1\begin{bmatrix}-\cos{t}+\sin{t}\\ 2\cos{t} \end{bmatrix}e^{2t}$$")

Figure 1: This is a simple mermaid graph.

Main classes of ODE models, ordered by increasing complexity

  • First order equations:

  • Michaelis-Menten equations

  • Auto-catalytic and Feedback Models

    • Positive activation
    • negative activation
  • Complex combination of several processes, pairing together multiple variables

Personalised medecine applied to Inflammatory Bowel Disease

1) Steady-state conditions

Note 1: Derive Steady-States Conditions of an ODE model

Computing the steady-state conditions for an ordinary differential equation (ODE) means finding the system’s equilibrium points, where the variables remain constant over time. Given a system of \(n=15\) ODEs (number of varying species in our model):

\[ \begin{cases} \frac{dx_1}{dt} &= f_1(x_1, x_2, \dots, x_{15}) \\ &\vdots \\ \frac{dx_{15}}{dt} &= f_n(x_1, x_2, \dots, x_{15}) \end{cases} \]

the steady-state conditions are obtained by solving system Equation 1:

\[ \begin{cases} f_1(x_1^*, x_2^*, \dots, x_{15}^*)& = 0 \\ & \vdots \\ f_n(x_1^*, x_2^*, \dots, x_{15}^*) &= 0. \end{cases} \qquad(1)\]

1) Steady-state implementation

healthy_model <- loadModel("../../models/team_2016_final_model_lo2016_2025_05_13.cps")
healthy_model_steady_state <- runSteadyState(
  model = healthy_model,
  calculate_jacobian = TRUE,
  perform_stability_analysis = TRUE,
  method = list("use_newton"=TRUE, accept_negative_concentrations = FALSE)
)
1
Run CoRC::loadModel() function to load the COPASI model in R.
2
We use the CoRC::runSteadyState function to run the steady-states…
3
… and for this task, several algorithms can be chosen, when no explicit solution can be derived directly. We used in our framework the well-known root-finding Newton-Raphson algorithm.

1) Steady-state results

Species

Value(g/cm^3)

M1 macrophage

1.17e-02

M2 macrophage

8.17e-03

Th1 cell

1.40e-01

Th2 cell

3.99e-03

Th17 cell

6.18e-03

Treg cell

2.98e-04

NA

2.01e-07

IL-2

1.07e-08

IL-4

3.36e-09

IL-21

7.77e-08

IL-6

1.46e-06

NA

2.72e-05

IL-10

3.01e-08

NA

1.64e-13

IL-12

3.64e-07

doi:10.1371/journal.pone.0165782.t004

Table 1: Steady-state concentrations of cytokines and immune cells, in healthy individuals in our model.

Original SS reported in Table 4, (Lo et al. 2016, 10).

Original SS reported in Table 4, (Lo et al. 2016, 10).

2) Sensitivity analyses


i) Sensitivity analysis Objectives

Sensitivity analysis aims to identify how changes in parameters affect model behaviour, identifying the most influential parameters on the outcome. Besides, these simulation studies contribute to prioritise where uncertainty matters most, and early capture bad representations of the model.


9 parameters were assumed in Lo et al. (2016) to have the strongest impact on the concentrations of T cell subsets. We retrieve their values, inferred from healthy individuals, using the CoRC::getParameters function (see Table 2):

Code
parameters_healthy <- getParameters(model = healthy_model,
                                  key = c("(Diff of M0 to M1).delta_m_cit",
                                          "(Diff of M1 to M2).sigma",
                                          "(Induction of T1 from M).s12",
                                          "(Proliferation of T1).s2",
                                          "(Induction of T2).s4",
                                          "(Induction of T17).s21",
                                          "(Induction of T17).s6",
                                          "(Induction of Tr).sb",
                                          "(Induction of Tr).s10")) |> 
  dplyr::select(-mapping)

flextable(parameters_healthy) |> 
  add_footer_row(values = "key is direct ID of the parameter in the model, reaction describes the biological mechanism associated with the value of the parameter, and value returns the constant value assumed for healthy individuals.", 
                 colwidths = 4) |> 
  bold(part = "header") 

key

name

reaction

value

(Diff of M0 to M1).delta_m_cit

delta_m_cit

Diff of M0 to M1

2.40

(Diff of M1 to M2).sigma

sigma

Diff of M1 to M2

24.00

(Induction of T1 from M).s12

s12

Induction of T1 from M

10.93

(Proliferation of T1).s2

s2

Proliferation of T1

1.23

(Induction of T2).s4

s4

Induction of T2

1.94

(Induction of T17).s21

s21

Induction of T17

156.17

(Induction of T17).s6

s6

Induction of T17

156.17

(Induction of Tr).sb

sb

Induction of Tr

14.02

(Induction of Tr).s10

s10

Induction of Tr

14.02

key is direct ID of the parameter in the model, reaction describes the biological mechanism associated with the value of the parameter, and value returns the constant value assumed for healthy individuals.

Table 2: The 9 reaction parameters evaluated for sensitivity analyses.

2) Sensitivity analyses


ii) Latin HyperCube Sampling, step I

Random Sampling for a \(d-\) random vector is performed as Equation 2:

\[ x_i^{(j)} \sim \mathcal{U}\left(x_{\min_j}^{(j)}, x_{\max_j}^{(j)}\right), \quad j = 1, 2, \dots, d, \quad i = 1, 2, \dots, N \qquad(2)\]

Figure 2: Random sampling, from Rustell (n.d.). With a low number of observations, the clustering pattern inherent to standard random sampling is clearly showcased.

Out of pure randomness, as shown in Figure 2, some regions of the bi-dimensional parameter space remain unsampled, while others are oversampled, as evidenced by visible clustering of points.

In Latin Hypercube Sampling (LHS):

  1. Sampling space is stratified into \(N\) equally-binned intervals for each dimension \(j\)
  2. Permutation is applied across all dimensions to avoid sampling correlation.
  3. Each sample is drawn from a distinct \(N-\) defined \(d-\)sized interval.
Figure 3: Latin hypercube sample, from Fig.3 Rustell (n.d.).

LHS provides a better Coverage of the Input Space and Improved Convergence, especially when the sampling size is large with respect to the number of draws, see Figure 3.

Main parameters:

  • \(x_i^{(j)}\) is the \(i\)-th sample drawn for the \(j\)-th dimension

  • \(d=9\) (number of parameters included in the sensitivity analysis)

  • \(N=5000\) the number of parameter distributions simulated

  • \(\mathcal{U}(x_{\min_j}^{(j)}, x_{\max_j}^{(j)})\) is the Uniform distribution, with bounds defined independently per dimension.

2) Sensitivity analyses


ii) Latin HyperCube Sampling, step II

n_samples <- 5000

parameters_sensitivity <- parameters_healthy$value
param_ranges <- tibble::tibble(parameters_name = parameters_healthy$key,
                               LB= parameters_sensitivity *0.8,
                               UB = parameters_sensitivity *1.2)


set.seed(20)
bootstrap_parameters <- lhs::randomLHS(n_samples, length(parameters_sensitivity))
for (i in seq_along(parameters_sensitivity)) {
  bootstrap_parameters[,i] <- qunif(bootstrap_parameters[,i],
                                    param_ranges$LB[i],
                                    param_ranges$UB[i])
}
apply(bootstrap_parameters,2, quantile)
1
Determine hyper-parameters of sensitivity analysis, such as the number of simulations \(N\), and lower and upper bounds within a \(\pm 20\%\) range.
2
We fix the random seed to ensure reproducibility of the sensitivity analyses.
3
The core function of the lhs package, lhs::randomLHS(), is only able to perform standard uniform sampling.
4
Apply the Inverse Transform Sampling Theorem (see Tip 1) to switch from a standard uniform distribution (\(U \sim \mathcal{U}(0, 1)\)) modelling to any bounded uniform distribution1
5
Apply the quantile function to verify the quantiles are rather close from their expected theoretical values from uniform sampling.

Tip 1: The Inverse Transform Sampling Theorem

The Inverse Transform Sampling Theorem allows generating any probability distribution from the standard uniform distribution, \(U \sim \mathcal{U}(0,1)\). Let \(X\) be a continuous random variable with cumulative distribution function (CDF) \(F(x)\) and \(F^{-1}\) its reciprocal (the quantile function), then:

\[ X = F^{-1}(U) \]

follows the same distribution as \(X\).


Specifically, by applying the following affine transformation \(X=a + (b - a) U\) on a standard uniformly distributed variable \(u\), \(X\) will follow this distribution: \(\mathcal{U}(a,b)\)

3) Partial Correlation Analyses


Standard Pearson correlation between two variables \(X\) and \(Y\) in Equation 3:

\[ \rho_{X,Y} = \frac{\text{Cov}(X, Y)}{\sigma_X \sigma_Y} \qquad(3)\]

  • \(\text{Cov}(X, Y) = \mathbb{E}[(X - \mu_X)(Y - \mu_Y)]\) is the empirical covariance matrix.
  • \(\sigma_X = \sqrt{\text{Var}(X)}\) and \(\sigma_Y = \sqrt{\text{Var}(Y)}\) are standard deviations.

Measures the strength of the linear relationship between \(X\) and \(Y\), without controlling for confounding variables.

The Partial Correlation measures the linear relationship between two variables, conditioned on the remaining set of variables \(\mathrm{Z}\), given by Equation 4:

\[ \rho_{X,Y \mid \mathbf{Z}} = -\frac{\Omega_{X,Y}}{\sqrt{\Omega_{X,X} \Omega_{Y,Y}}} \qquad(4)\]

  • \(\mathbf{\Omega} = \mathbf{\Sigma}^{-1}\) is the precision matrix.
  • Non Null off-diagonal terms of the precision matrix reflect direct, causal connections between two variables.

Partial correlation enables more straightforward interpretation by revealing truly causal relationships, while discarding spurious connections.

Practical use case: in Figure 4,

  • The standard Pearson correlation would have certainly observed a significant correlation between variables \(A\) and \(C\).
  • This spurious connection between \(A\) and \(C\) due to the confusing effect of \(B\) is clipped using partial correlation instead.
Code
digraph G {
    layout=neato
    A -> B;
    B -> C;
}
G A AB BA->B C CB->C
Figure 4: Simple example of chain rule associations.

4) Rank Transformation and Scatter Plots


Original Fig.2. reproduced from (Lo et al. 2016, 11)

Fig.2. resulting from our own sensitivity analyses.
Figure 5: Scatter plots of rank-transformed concentrations of \(T_1\) and \(T_2\) immune cells, against 3 variations of parameters. We can depict at least 2 errors or biological inconsistencies from original simulations: PRCC is not strictly bounded between -1 and 1 (while it should have been normalised by the variances), and from (zenewicz2009timm?), and ODE Equations 12 and 13 from Lo et al. (2016), we would have expected a Th2 inhibition by Th1, in other words, a negative PRCC coefficient for factor \(\sigma_{12}\).

BioModels: a repository of biological models

BioModels: Metadata description

Around half of the publicly avalaible models are tagged manually-curated, a ratio slightly increasing over time.

Main reasons underlying failure of reproducing models.
Figure 6: BioModels key features.

BioModels submission flowchart

From submission to publication of a model in BioModels. Reproduced from Fig4, from Malik-Sheriff et al. (2020)

From raw SBML format to curated status. Reproduced from Fig2.A, from Malik-Sheriff et al. (2020).
Figure 7

Conclusion

Take-Aways:

  • BioModels, a data platform for storing curated biological models.
  • However, no automated running of code snippets, nor mandatory access to original data
  • Hard-coded tools, such as COPASI, enable quick analyses for non-expert audience
  • But they’re limited in terms of functionalities and flexibility, to be used with caution.

This is a Quarto RevealJS Presentation. To learn more, visit https://quarto.org/docs/presentations/revealjs.


Bibliographic References

Baker, Monya. 2016. ‘1,500 Scientists Lift the Lid on Reproducibility’. Nature 533 (7604): 452–54. https://doi.org/10.1038/533452a.
Begley, C. Glenn, and Lee M. Ellis. 2012. ‘Raise Standards for Preclinical Cancer Research’. Nature 483 (7391): 531–33. https://doi.org/10.1038/483531a.
Cokelaer, Thomas, Sarah Cohen-Boulakia, and Frédéric Lemoine. 2023. ‘Reprohackathons: Promoting Reproducibility in Bioinformatics Through Training’. Bioinformatics 39 (June): i11–20. https://doi.org/10.1093/bioinformatics/btad227.
Initiative, The Brazilian Reproducibility, Olavo Bohrer Amaral, Clarissa Franca Dias Carneiro, Kleber Neves, Ana Paula Wasilewska Sampaio, Bruna Valerio Gomes, Mariana Boechat de Abreu, et al. 2025. ‘Estimating the Replicability of Brazilian Biomedical Science’. 25 April 2025. https://doi.org/10.1101/2025.04.02.645026.
Lo, Wing-Cheong, Violeta Arsenescu, Razvan I. Arsenescu, and Avner Friedman. 2016. ‘Inflammatory Bowel Disease: How Effective Is TNF-α Suppression?’ PloS One 11 (11): e0165782. https://doi.org/10.1371/journal.pone.0165782.
Malik-Sheriff, Rahuman S, Mihai Glont, Tung V N Nguyen, Krishna Tiwari, Matthew G Roberts, Ashley Xavier, Manh T Vu, et al. 2020. BioModels—15 Years of Sharing Computational Models in Life Science’. Nucleic Acids Research 48 (D1): D407–15. https://doi.org/10.1093/nar/gkz1055.
Open Science Collaboration. 2015. ‘Estimating the Reproducibility of Psychological Science’. Science 349 (6251): aac4716. https://doi.org/10.1126/science.aac4716.
Rustell, Michael. n.d. ‘Knowledge Extraction and the Development of a Decision Support System for the Conceptual Design of Liquefied Natural Gas Terminals Under Risk and Uncertainty.’ University of Surrey. Accessed 14 May 2025. https://openresearch.surrey.ac.uk/esploro/outputs/doctoral/Knowledge-extraction-and-the-development-of/99511011302346.

About

License rights

© 2025 Bastien CHASSAGNOL

Opinions expressed are solely my own and do not express the views of the researchers I am associated with.

This work is released under the MIT license.

Made by Quarto Presentations RevealJs